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Abstract The focus of this work is on an alternative implementation of the iterative ensemble smoother 
(iES). We show that iteration formulae similar to those used in [2E can be derived by adopting a 
regularized Levenberg-Marquardt (RLM) algorithm ||T4] to approximately solve a minimum-average-cost 
(MAC) problem. This not only leads to an alternative theoretical tool in understanding and analyzing the 
behaviour of the aforementioned iES, but also provides insights and guidelines for further developments 
of the smoothing algorithms. For illustration, we compare the performance of an implementation of the 
RLM-MAC algorithm to that of the approximate iES used in [3] in three numerical examples: an initial 
condition estimation problem in a strongly nonlinear system, a facies estimation problem in a 2D reservoir 
and the history matching problem in the Brugge field case. In these three specific cases, the RLM-MAC 
algorithm exhibits comparable or even better performance, especially in the strongly nonlinear system. 


Introduction 

For data assimilation problems there are different ways in utilizing the available observations. While 
certain data assimilation algorithms, for instance, the ensemble Kalman filter (EnKF, see, for example, 
EE), assimilate the observations sequentially in time, other data assimilation algorithms may instead 
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collect the observations at different time instants and assimilate them simultaneously. Examples in this 
aspect include the ensemble smoother (ES, see, for example, 0]) and its iterative variants. 

The EnKF has been widely used for reservoir data assimilation (history matching) problems since 
its introduction to the community of petroleum engineering [261 . The applications of the ES to reservoir 
data assimilation problems are also investigated recently (see, for example, [32]). Compared to the EnKF, 
the ES has certain technical advantages, including, for instance, avoiding the restarts associated with 
each update step in the EnKF for certain reservoir simulators (e.g., ECLIPSE© [1]) and also having 
fewer variables to update. The formal benefit (avoiding restarts) may result in a significant reduction of 
simulation time in certain circumstances [33], while the latter (having fewer variables) can reduce the 
amount of computer memory in use. 

To further improve the performance of the ES, some iterative ensemble smoothers (iES) are suggested 
in the literature, in which the iterations are carried out in the forms of certain iterative optimization 
algorithms, e.g., the Gaussian-Newton mm or the Levenberg-Marquardt method 
or in the context of adaptive Gaussian mixture (AGM, see [21] )• In [131122] . the iteration formulae are 
adopted following the regularized Levenberg-Marquardt (RLM) algorithm in the deterministic inverse 
problem theory (see, for example, EES])- Essentially these formulae aim to find a single solution of 
the inverse problem, and the gradient involved in the iteration is obtained either through the adjoint 
model [13], or through a stochastic approximation method [22]. While in the iteration formula is 
derived based on the idea that, the final result of the iES should be equal to the estimate of the EnKF, 
at least for linear systems with Gaussian model and observation errors. Consequently, this algorithm is 
called ensemble smoother with multiple data assimilation (ES-MDA for short). On the other hand, in 
[2] an iteration formula is obtained based on the standard Levenberg-Marquardt (LM) algorithm. By 
discarding a model term in the standard LM algorithm, an approximate iteration formula is derived, 
which is similar to that in [Bj. For distinction, we call it approximate Levenberg-Marquardt ensemble 
randomized maximum likelihood (aLM-EnRML) algorithm. 


1 These optimization algorithms are also applied to iterative EnKFs, e.g., in some of the aforementioned works. 
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In this work we show that an iteration formula similar to those used in the ES-MDA and RLM-MAC 
can be derived by applying the RLM to find an ensemble of solutions via solving a minimum-average-cost 
problem m- The gradient involved in the iteration is obtained in a way similar to the computation 
of the Kalman gain matrix in the EnKF. This derivation not only leads to an alternative theoretical 
tool in understanding and analyzing the behaviour of the aforementioned iES, but also provides insights 
and guidelines for further developments of the iES algorithnjf|. As an example, we derive an alternative 
implementation of the iterative ES based on the RLM algorithm. Three numerical examples are then 
used to illustrate the performance of this new algorithm and compare it to the aLM-EnRML. 


Methodologies 


ES-MDA and aLM-EnRML 


Let m denote an m-dimensional reservoir model that contains the petrophysical properties to be esti¬ 
mated, g the reservoir simulator, and d° a p-dimensional vector that contains all available observations 
in a certain time interval. Here d° is assumed to contain certain measurement errors, with zero mean 
and covariance C^. In the context of iterative ensemble smoothing, suppose that there is an ensemble 
M‘ = of N e reservoir models available at the ith iteration step, then an iES updates M’ to its 

counterpart M !+1 = {ml +1 }^f 1 at the next iteration step via a certain iteration formula, which is the 
focus of our discussion below. For convenience of discussion later, let us define the following square root 
matrix S* ra with respect to the model m (model square root for short): 


CJl - 


VNe - 1 


N e 


N e 


i=i 


(1) 


in the sense that the product S l m (S^) 7 is equal to the sample covariance of ML 

The idea in the ES-MDA is to make the final iES estimate equal the estimate of the EnKF, at least for 
linear systems with Gaussian model and observation errors (some thoughts are also provided in Appendix 
El in situations with nonlinearity and/or non-Gaussianity). To this end, an iteration formula is derived 

2 For instance, instead of using the RLM algorithm, one could apply other types of regularized inversion methods, e.g., 
the regularized Gaussian-Newton or conjugate gradient algorithm 1711151 . to solve the minimum-average-cost problem. This 
will thus lead to more types of iterative ES algorithms. 
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in terms of [(] 


" i+1 = mj + Si (S i d ) T (s j (S > d ) T + Y C d ) (d° - g (mj)) , i = 1, 2, ■ ■ ■ , J; 


VN e - 1 L 

JHl-lEsK) 


g(mi)-g(mt),... ,g(m^ e )-g(mj) 


J = 1 


The total number of iteration steps / is chosen before the iteration starts; d° (j = 1, 
perturbations of d°, similar to those in the EnKF; and Y (* = !,•■■ , I) satisfy 


i 1 

7 ® > 1 and —7=1. 


( 2 ) 

(3) 

(4) 

■ • , N e ) are the 

(5) 


For convenience, we also call data square root hereafter, in the sense that S l d (S l d ) T is equal to the 

sample covariance matrix of the simulated data {g (m*)}^ 1 . 

In the Levenberg-Marquardt ensemble randomized maximum likelihood (LM-EnRML) method [3], 
history matching is recast as a weighted least-squares problem, in the sense that each model, say mj (the 
j th ensemble model at the ith iteration step), is updated to the one mj +1 at the next iteration step by 
solving the following minimization problem: 


argmin (d° — g (m)) T C d 1 (d° — g ( m )) + ( m — mj) T ( < -‘m) 1 (m - mj) , (6) 

m 

with C° m = S° m (S^) T , and given by Eq. (HD. Applying the Levenberg-Marquardt algorithm [23] to 
©, one has the following iteration formula to update the model [3j 


m 


-Am) 


(7) 


4mj = [(C°,) 1 + (Gj) T C7 1 G5 


] 1 [(G;-) T CJ 1 (dJ-g(mj))-(C° l ) ‘(m’-mj)] 


[(l + V)( C j„) 1 + (G}) T C^ 1 Gj] 1 [(Gj) T c7(d°-g(mj))-(C^) 1 (mj - m°)] 


r p 

where C z m = (Sj„) , Gj is the Jacobian of g evaluated at mj, and A* a positive scalar. In the course 

of deriving the third line of Eq. 0, the first (C^) 1 is approximated by (1 + A 1 ) (C^) 1 , in accordance 
to the Levenberg-Marquardt algorithm [3]. In large-scale problems, it can be expensive to evaluate the 
term (CjjJ -1 (mj - mj), because the model size m is usually much larger than the observation size p 
(i.e., to p ). Hence discarding the term (CjjJ 1 (mj — mj), using the Sherman-Morrison-Woodbury 
formula [32] and some algebra, one has 


A+i 


+ CjjGjf [GjCjjGjf + (1 + V) C d ] 1 (d° - g (mj)) . 


( 8 ) 


As to be shown later (see Eq. (ED), with a suitable choice for A*, the iteration formula in Eq. 0 is 
consistent with the one derived from the RLM algorithm. Under suitable technical conditions (see, for 
example, DU), mj will converge to a solution of the equation g (m) = dj as * —► + 00 . 










Title Suppressed Due to Excessive Length 


5 


Following the EnKF, one may introduce a further approximation to the gain matrix in Eq. ©, such 


that 


(si(Sj) T + (l + A l )C d ) 


Cm(GJ) J [GJCJJGJ)' + (1 + A 1 ) Cd 
Therefore, with the above approximations, one has the following final iteration formula 

= mj + (S^) T (S* (S^) T + (1 + V) Cd) _1 (d° - g (mj)) . 




(9) 


( 10 ) 


Comparing Eqs. (H and ©> it is clear that the iteration formulae become identical when (1 + A*) in 
Eq. m ^ replaced by 7 * in Eq. ©. 

For distinction, we call the algorithm in Eq. m approximate Levenberg-Marquardt ensemble ran¬ 
domized maximum likelihood (aLM-EnRML) method. In the aLM-EnRML, starting from an initial value 
A 0 , the subsequent values of A* will decrease if the average data mismatch with respect to the models is 
reduced, otherwise the values of A* will increase instead. The iteration process will stop if the maximum 
iteration number is reached, or if the relative change of the average data mismatch in two consecutive 
iteration steps is lower than a given threshold. 


The regularized Levenberg-Marquardt (RLM) algorithm and an extension 

In the sequel we first introduce the regularized Levenberg-Marquardt (RLM, as used in [T31I22] ) algorithm 
in the context of deterministic inverse problems theory, and then consider an extension of the RLM 
algorithm as an ensemble data assimilation method. In the course of deduction, the similarities between 
the ES-MDA/aLM-EnRML and the extended method in the current study will become clear, while some 
implications in the extended method will also be discussed. Note that in a recent work na, the author also 
considers using the RLM algorithm in the context of ensemble smoothing. Compared to m , the current 
study has different focuses, e.g., in terms of the way in approximating the gradient of the cost function of 
the minimization problem, and the parameter rule used in the RLM algorithm (see (1161) later Jf| Readers 
are referred to [ 12 ] for more detail. 

3 To summarize, mmm and the current work are all related to the RLM algorithm. mm use the parameter rule 


based on, for example, the work m, while m and the current work use the parameter rule based on UU- Although 
different in the concrete forms, both parameter rules are proven to lead to local convergence of the RLM algorithm mm- 
In addition, I13II22| both aim to find a single solution, but employ different methods in approximating the gradient of the 
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In the conventional deterministic inverse problem theory (see, for example, [71 1151 1. one aims to find 
a single set of parameters (e.g., a reservoir model m) whose simulated output (e.g., g(m)) matches the 
observation d°. Intuitively, one can search for such a model by solving the following weighted least-squares 
problem 

argmin (d° - g (m)) T C ^ 1 (d° - g (m)) , ( 11 ) 

m 

which aims to minimize the data mismatch (d° — g (m)) T C ^ 1 (d° — g (m)), where C ^ 1 is the weight 
matrix associated with the data mismatch term. For convenience of discussion later, let us define ||d ||2 = 
Vd T d and ||d||c d = \Jd r C~^ d for a vector d, then it is clear that (d° — g ( m )) Tc d 1 ( d °-g( m )) = 
||d°-g(m)||^ = ||C; 1 / 2 (d°-g(m))|||. 

When the dimension of the observation space is lower than that of the model space (as is often the 
case in practice), (1111) is an under-determined and ill-posed problem, which has some undesirable features, 
e.g., non-uniqueness of the solution and potentially large sensitivity of the solution to the observation, in 
the sense that even a very small change in the observation might lead to a large change in the solution [7]. 
To mitigate the above problems, in deterministic inverse problems theory it is customary to introduce a 
certain regularization term to ED- Here we consider Tikhonov regularization, in which a regularization 
term, in terms of (m — m fc ) T C ” 1 (m — m fc ), is introduced to ED- such that the weighted least-squares 
problem becomes 

argmin (d° - g (m)) T Cj 1 (d° - g (m)) + 7 (m - m 6 ) C ” 1 (m - m 6 ) , (12) 

where 7 is a positive scalar, m b denotes a background model, and C ” 1 is the associated weight matrix. 
The choice of 7 affects the relative weights assigned to the data mismatch term and the regularization 
term, and thus influences the solution of ED- For instance, if 7 is relatively large (e.g., tending to 00 
in the extreme case), then the obtained solution will approach m b . On the other hand, if 7 is relatively 
small (e.g., tending to 0 ), then (fill) and (fTTl) tend to coincide with each other, and the solution of 
ED can be considered as an approximation to one of the solutions of (1111) . Apart from 7 , the matrices 
C m and Cd would also affect the relative weights between the regularization and data mismatch terms. 
If, in addition, C m and are chosen to be the error covariance matrices of the background and the 

cost function, as discussed previously. In contrast, m and the current study target for multiple solutions in the context of 
iES, and differentiate each other in terms of the ways in constructing the data square root matrices (e.g., m uses the same 
way as in the EnKF). In Appendix [X] we also provide an alternative point of view to interpret the iES algorithm derived 


in the current study. 
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observations, respectively, then the correlations of the variables in the model and observation spaces, 
respectively, are also incorporated in the minimization problem. 

In general, one may wish to choose a 7 such that the data mismatch of the resulting solution is 
comparable to the noise level of the observations mm- The true noise level is often unknown, therefore 
in practice it may be replaced by, for instance, a threshold value proportional to p (i.e., the size of d°). 

Readers are referred to, for example, mmm for the rationales behind this choice. 

For linear systems, there are straightforward ways for one to construct a solution that has the desired 
data mismatch (see, for example, [23] ). For nonlinear systems, however, one often has to rely on a certain 
iterative algorithm to find the solution. One of such methods, called the regularized Levenberg-Marquardt 
algorithm following |14] (or the regularizing Levenberg-Marquardt algorithm following [IT]), constructs 
a sequence of model {m 1 } (i = 1 , 2 , ■ • ■) by solving a linearized weighted least-squares problem 

argmin(d° - g (m 4 ) - G;(m i+1 - nT)) 7, C^ 1 (d° - g (m 4 ) - Gi(m i+1 - m ,: )) (13) 

m i + 1 

+7 i (m ,:+1 - m i+1 ’ i ’) T C" 1 (m 4+1 - m i+1 ’ b ) , 


at each iteration step. In (fT3l) . G; is the Jacobian of g evaluated at m', such that g (m*) +Gi(m I+1 -m') 
represents a first order Taylor approximation to g (m ,+1 j; 7 * is a positive scalar as in m, but now 
changes over the iteration; m 4+1 ’ b is the background model that may also be adaptive with the iteration. 
A convenient choice can be m I+1,b = m 1 , meaning that we want the new model m I+1 not to be too far 
away from the previous model, such that the approximation g (m' +1 ) ss g (m ! ) + Gi(m !+1 — m l ) can 
be roughly valid [22]. In addition, the choice m i+1,b = m' also simplifies the iteration formula, as will be 
shown below. 

The solution of the weighted least-squares problem (THU) is given by 


m i+1 =m i + ZIm i (14) 

4W = [VC- 1 + (G*) T C^G 4 ] -i [(G’f C^ 1 (d° - g (m 4 )) - VC- 1 (m i - , 

which is similar to Eq. 0 (e.g, by replacing 1 + A* therein with 7 *). By letting m i+1,b = m' and with 
some algebra, one has the simplified iteration formula 


- i+1 =m 1 + C m (G i ) T |G i C m (G l ) T + 7 i C d | ‘ (d° - g (m')) 


(15) 


which is thus similar to Eq. (in. One may also let the weight matrix C m be adaptive with the iteration, 
which will be considered later in the context of iES. Intuitively, provided that the step size Am 1 is small 
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enough (which can be controlled by 7 *) such that the first order Taylor expansion of g(m i+1 ) around m 1 
is a roughly valid approximation of g(m !+1 ), one has 


Hgfm*"*" 1 ) - d°|fe « ||g (m 4 ) - d° + G* (m 1+1 - in 1 ) ||®. 


= ||[I — G'C m (G') T G i C m (G i ) T + 7 i C d ] (g( mi ) — d°) 


= || 7 l C d |G*C m (GV + 7 1 C d J (g(m‘) - d°) \\ 2 Cd 

= (C d 1/2 (g(m*) - d°)) T [(CJ^G 4 ) C m (CJ 1 / 2 G i ) T + 7 * Ip] * (c d 1/2 (g(m 4 ) - d°)) 

= H7 1 [(cj^G 4 ) C ro (CJ 1 /a G 4 ) T + 7 “ I p 
i-i / 2 rJT . 


,-1/2 
J d 

< l|g( m *) — d°||Q , 


C d 1/2 (g(m 4 )-d»)) 


< ||7* [(CJ : ^G ') C ra (C- 1 "G I f + 7 l I P | ||£ l|g( m “) - d 0 ||^ rf 


since ||7 Z C d 1/2 G 4 ) C m (C d 1/2 G i ) T + 7 l IpJ ||| < 1, where I p denotes the p-dimensional identity matrix, and || A||2 

represents the spectral norm of the matrix A. This implies that the data mismatch ||g(nT~* _1 ) — d°||^ d at the 

(i + l)-th iteration tends to be no larger than the one ||g(m 2 )- d °llc d at the previous iteration. 

In terms of the choice of 7 % the following parameter rule from e 

7 ° > 0, 

7 l+1 = p7*, with 1 /r < p < 1 for some scalar r > 1, ( 16 ) 

lim 7 l = 0 , 


can be used, in which the scalar sequence { 7 *} gradually reduces to zero as i tends to + 00 , while the 
presence of the lower bound 1 jr for the coefficient p aims to prevent any abrupt drop-down of 7 * to zero. 
When Eq. (flbl) is used in conjunction with (flfill . it can be analytically shown that the data mismatch of 
the sequence of models {m 1 } converges to zero locally as i —> +00 and the observation noise level tends 
to zero, provided that the equation g(m) = d° is solvable and some other conditions are satisfied E- 
Instead of using OSD, one can adopt other parameter rules that also lead to local convergence of the data 
mismatch under similar technical conditions (see, for example, lllj ). 

When the observation noise is present, however, it may not be desirable to have a too small data 
mismatch in order to prevent over-fitting the observations. In general, one may wish to let the iteration 
stop when the data mismatch of the current iterated model is comparable to the noise level, which is 
often referred to as the discrepancy principle mm- Of course, in practice the true noise level is typically 
unknown. As a result, one may let the iteration process Eq. (1151) stop when either of the following three 
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conditions are satisfied: (a) the data mismatch of the iterated model is lower than a pre-set threshold 
f3^P f° r the first time for a given positive scalar /3 U , in light of the discussions in [rmresiresiisni : or (b) the 
maximum iteration number is reached; or (c) the relative change of the data mismatch in two consecutive 
iteration steps is lower than a given threshold (0.01% in our implementation). Note that, with the stopping 
condition (a), the total number of (actual) iteration steps can be less than the maximum iteration number, 
for instance, when the RLM reaches the threshold f3^p earlier than up to the maximum iteration number. 
Conditions (b) and (c) are the same as those in the aLM-EnRML, and are mainly introduced to control 
the runtime of the iteration process. These settings will be applied to the aLM-EnRML and our extension 
scheme in the experiments later, in order to make the comparison of both algorithms be conducted under 
the same experimental settings as far as possible. 

In the practical implementation of the RLM algorithm (and our extension scheme below), there are a 
few issues that need to be taken into account. Firstly, in the theoretical analysis of the RLM algorithm (see, 
for example, Emailing), it typically assumes that the initial model stays sufficiently close to a solution of 
the equation g(m) = d°, which is not necessarily true in reality. In addition, as can be seen in our previous 
discussion, the validity of the first order Taylor approximation g (m‘ +1 ) ss g (m*) + Gi(m‘ +1 — m 1 ) is 
essential to the derivation of the RLM algorithm. In practice, it is often not a trivial issue to obtain the 
Jacobian G^. If the exact Jacobian G, is not available, then one must evaluate it based on a certain 
approximation scheme, which often has to struggle between the computational efficiency and accuracy. 
Putting these factors together, in practice it is likely that the iteration process, Eq. m may occasionally 
lead to higher data mismatch for the current iterated model, in comparison to that of the previous one, 
due to, for example, a too large step size that invalidates the first order Taylor approximation. In such 
cases, following [3], a simple remedy can then be to increase the parameter y 1 to some extent and re- 
execute the iteration in Eq. (11511 (which will then have a smaller step size). This strategy is essentially 
similar to the back-tracking line search method in optimization theory (see, for example, 128 ] ) and works 
well in our numerical studies. 

In the remainder of this section we consider an extended RLM scheme in the context of ensemble 
smoothing. As in the aLM-EnRML, instead of examining the data mismatch of a single model, we are 
more interested in studying how the average data mismatch of the ensemble models changes over the 
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iteration step. To this end, the weighted least-squares problem in (fT!?l) is replaced with the following 
minimum-average-cost (MAC) problem, namely, 


. a ^r e w e £ [<5 - g H +1 )) T c * 1 ( d ° - g B +1 )) 

{"V b=l j=l 


+T'(m‘ +1 -m; 


m } +1 - m 7 


(17) 


where C( n = (S^j (S^ will be specified later), and d", m* and mj +1 are the same as those 
defined previously. For distinction, we call the resulting iES regularized Levenberg-Marquardt algorithm 
based on minimum-average-cost (RLM-MAC for short). Appendix [Al also includes some thoughts that 


aim to interpret the RLM-MAC from the point of view of an expectation-maximization algorithm [5]. 

To solve 03. it is necessary to linearize the forward model g as in the RLM algorithm. Note that in 
mi there are N e simulated observation terms g (m* +1 ) (j = 1,2,--- , N e ). To reduce the computational 
cost in evaluating the Jacobian matrices, one may choose to linearize the forward model g, for each 
ensemble member m* +l , around a “common” point, say m*, such that 


g 



; K) + g* 



(18) 


where G). is the Jacobian of g at m*. Inserting Eq. (TTH1) into (fI71) and with some algebra (see Appendix 
03, one has the following iteration formula 


m* +1 = raj + 4m) 


zlmj = S i rn (G i c S i m ) T [(G(sjj(G(S(J T + 7 ‘ C d ] (d° - g (mj) - G l c (mj - m*)) , 
for j = 1, • • • , N e . If one lets 


C* - 


1 


VN e - 1 


m, — m,. 


, m„ - m„ 


and adopts the approximation G), (m) — mj) « g (mj) — g (mj), then Eq. (fTHl) reduces to 
m) +1 = mj + SjJSjf (Si(SSr + 7* G.y 1 (d° - g (mj)) , 

where 


Si, = 


1 


-1 


[g (ml) - g (mj) , • ■ ■ , g (m l N ) - g (mj)] . 


(19) 


( 20 ) 


( 21 ) 


( 22 ) 


Compared to the iteration formulae of ES-MDA and aLM-EnRML (e.g., Eq. (USD), the iteration formula 
of RLM-MAC in general differs in the way of constructing the square root matrices Sj„ and S l d . Possible 
choices for mj could be, for instance, the mean m' or a member (e.g., the one closest to m‘ ) of the 
ensemble {m*}j^ 1 . In the current work we let mj = m 1 . 
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Under the choice mj, = m 1 , coincides with S ? m , therefore the resulting RLM-MAC algorithm 
mainly differs from the aLM-EnRML in the data square root matrix. In general, and may have 
different “centering points”: in Eq. m, the simulated observations g (m*) in center around g (m*), 
while in they center around g (in') instead. In certain circumstances, e.g., when the forward model 
g is linear or weakly nonlinear, g (m*) and g (m‘) may be close to each othej^. In some other cases, 
the distinction between g (m') and g (m l ) may be more substantial, and the behaviour of aLM-EnRML 
and RLM-MAC may thus become more diverged, as will be shown later. In terms of the performance 
comparison between the aLM-EnRML and RLM-MAC, however, we stress that, although one may see 
the relative superiority of one algorithm over the other in some specific cases, it is expected that the 
conclusion might not be valid in a broader sense, following the “no free lunch theorems for optimization” 
m- Intuitively, since both the aLM-EnRML and RLM-MAC are local optimization algorithms and they 
may have different search directions, these two algorithms might end up with different local optima in 
general, and the relative superiority between these two algorithms may thus vary from case to case. 

Finally we note that, in the course of solving (El). we have made use of the first order Taylor ap¬ 
proximation around a common point m* in order to reduce the computational cost. In this regard, a 
theoretically more sound - yet computationally more expensive - strategy would be to linearize the for¬ 
ward model g, for each ensemble model, around a certain point in its neighbourhood. A trade-off strategy 
may also be employed, e.g., by clustering the ensemble models into a small number of groups, and then 
linearizing g around a common point for each group. In addition, in the context of solving a MAC prob¬ 
lem, it is also possible for one to explicitly incorporate into ED additional terms that impose certain 
constraints onto the statistics of the ensemble members or the corresponding simulated observations, 
such that the iteration step may be in a form beyond the Kalman update formula as in Eq. ED- An 
investigation in this aspect will be carried out in the future. 


4 


In particular, when g is linear, g 



g (m l ), such that the RLM-MAC and aLM-EnRML become identical. 
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Numerical examples 


After discussing the similarities and differences between the ES-MDA, aLM-EnRML and RLM-MAC in 
the previous section, we focus on demonstrating potentially different behaviour of the ensemble smoothers 


that is resulted from using different “centering points” (g (m® ) or g (m * 1 )) to construct the data square 
root matrix. 

As noted previously, in the ES-MDA, the total number I of iteration steps is decided before the 
iteration process starts, and the parameters 7 ® are dependent on /, in light of the constraint (0. The 
aLM-EnRML and RLM-MAC do not have the above features in general: Under the stopping conditions 
(a-c) in the previous section, the total numbers of the iteration steps of the aLM-EnRML and RLM-MAC 
depend on the threshold /3and the convergence speed of the algorithms in specific problems, subject 
to the constraint of maximum iteration number. In addition, the parameters 7 ® are chosen according to 
a different rule, which is independent of I or the maximum iteration number in general. Therefore, for 
our purpose here, we confine ourselves to the comparison between the aLM-EnRML and the RLM-MAC. 
In the implementations, both algorithms apply the truncated singular value decomposition (TSVD) to 


and retain the leading 


the data square root matrices (normalized by C d 1 ^ 2 ) in their iteration formulae®, 
eigenvalues which sum up to at least 99% of the “total energies” (see [51155] for the details). In the 
experiments below we first examine the performance of both algorithms in a strongly nonlinear system, 
and then apply both algorithms to a reservoir facies estimation problem and the Brugge field case. 


An initial state estimation problem in a strongly nonlinear system 

Here we adopt the strongly nonlinear 40-dimensional Lorenz 96 (L96) model [TS] to demonstrate the 
potentially different behaviour between the aLM-EnRML and RLM-MAC. The governing equations of 

5 In doing so, the iteration formulae, e.g., that in Eq. (1211) . are re-expressed as m® 1 1 = m® + 

1 ^ 2 SJ i ) T ((C d 1/2 S^(C d 1/ "S^) T + 7 1 I p ^ C d 1 ^ 2 (A° — g (Ai®.^ , where l p is the p-dimensional identity matrix. 

The TSVD is applied to the product ( C d 1 2 S®()(C f/ 1 / 2 which can have certain numerical benefits, e.g., mitigating 
the issue of different orders of amplitudes in the measurement data. 
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the L96 model are given by 

dxb 

-jj- = (Xk+i - Xk- 2 ) Xk -1 - X k + F, k = 1, • • • , 40. (23) 

I 11 Eq. (E3l) the variables Xk are defined in a cyclic way such that x_i = X 39 , Xq = X 40 and X 41 = X\. 
The L96 model is designed to mimic baroclinic turbulence in the midlatitude atmosphere. Because of 
its strongly nonlinear nature, this model is often employed to test the performance of data assimilation 
algorithms in the weather data assimilation community. 

In this work we let F = 8 , with which the L96 model is chaotic. The L96 model is discretized 
through the fourth-order Runge-Kutta method with a constant integration step of 0.05 (corresponding to 
6 hours in real time, see nn). To initialize the model, we draw an initial state vector at random from the 
multivariate normal distribution ./V(0, 140), where I40 is the 40-dimensional identity matrix. To avoid the 
transit effect, the model states within the first 500 integration steps are discarded, and data assimilation 
starts after this “spin-up” period. For reference, hereafter we re-label the model state at time step 500 
as the initial state, which has an impact on subsequent model states and their observations. In the 
experiment the initial state vector is assumed unknown, and will be estimated based on the observations 
at subsequent time steps. 

In the experiment the assimilation time window is 10 days, and the observations are available every 
24 hours. At each observation time instant (e.g., 24h, 48h etc.), the synthetic observations are generated 
by first applying the nonlinear function f(x) = x 3 /5 to the odd number state variables X\,X 3 , ■ ■ ■ , 2:39 
(overall 20 variables), and then adding a sample from the normal distribution N( 0,1) to each function 
value f(xk ) (k = 1,3, • • • , 39). The total number of observation elements in the assimilation time window 
is thus 10 x 20 = 200 . 

The initial background ensemble is drawn from a multivariate normal distribution whose mean and 
covariance are obtained from the “climatological” statistics of the L96 model in a long free-run with 
100,000 integration steps (see, for example, [22] for the details). The ensemble size of the aLM-EnRML 
is 100. Note that in the RLM-MAC one also needs to compute g (m‘), the simulated observation of the 
ensemble mean, apart from those of the ensemble members. Because of this, we let the ensemble size of 
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the RLM-MAC be 99, such that the aLM-EnRML and RLM-MAC have the same number of simulated 
observations at each iteration step. The maximum number of iterations is 100, and the coefficient /3 U = 2. 

In light of the requirement that the first order Taylor approximation be roughly valid, it is suggested 
to start with a relatively small step size at the beginning [ 22 ]- As a result, the initial value 7 0 should 
be relatively large (but not too large in order to reduce the total number of iteration steps). In our 
implementation, we let 7 * = a 1 y / trace(S^(S^) T )/A e , where trace(*) calculates the trace of a matrix, 
a° = 1 and a l+1 = 0.9 a* if the average data mismatch of the ensemble is reduced, otherwise the ith 
iteration step is repeated, with a 1 = 2 a 1 , for maximum 5 times, following [2- In the subsequent two other 
case studies, the same way is applied to tune a* adaptively. We stress that this should not be considered 
as the optimal rule. Instead, it is expected that one may improve it for specific applications in general. 

As the dimension of the L96 model is relatively low, we can afford to repeat the experiment 100 times. 
In each repetition, the initial background ensemble and the associated observations are drawn at random. 
In the same repetition run, the aLM-EnRML and RLM-MAC have the same initial background ensemble, 
the associated observations and the parameter rule in choosing 7 *. 

Figures []]-[3] report the data mismatch and RMSEs of the aLM-EnRML and RLM-MAC, obtained 
from 100 repetitions of the experiment|f|. In terms of the data mismatch, from Figures |T]and HI one can see 
that for the aLM-EnRML, its range is roughly between 10 3 and 10 6 . Overall, 2% of the data mismatch 
lies in the interval [10 3 ,10 4 ], 51% in the interval [10 4 ,10 5 ], and 47% in the interval [10 5 ,10 6 ]. For the 
RLM-MAC, the range of the data mismatch is roughly between 10 2 and 10 6 . there are 11% of the data 
mismatch contained in the interval [10 2 ,10 3 ], 13% in the interval [10 3 ,10 4 ], 51% in the interval [10 4 ,10 5 ], 
and the remaining 25% in the interval [10 5 ,10 6 ]. This suggests that in some cases the RLM-MAC will have 
data mismatch lower than the minimum one of the aLM-EnRML. Similar results can also be observed 
in terms of the RMSEs (Figures [T] and [3]). As can be seen in Figure (TJ for the aLM-EnRML, its RMSEs 
seem to follow a uni-modal distribution, with its (single) peak around 3. On the other hand, for the RLM- 
MAC, the distribution of its RMSEs tends to be more flat, with one of its peaks being closer to zero. By 
counting the numbers of the RMSEs in the sub-intervals (e.g., [0,1], [1, 2] etc) of the range [0, 6 ], the pie 

6 Here the data mismatch is calculated according to GIli, and is averaged over all ensemble members. The RMSE of a 
m-dimensional estimate v with respect to its truth w tr is give by ||v — v tr \\ 2 / y/rn, where || • ||2 denotes the ^ 2 -norm. 
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charts in Figure [3] indicate that, the percentages of the RMSEs of the aLM-EnRML in the sub-intervals 
(from low to high values) are 2%, 18%, 34%, 31%, 12% and 3%, respectively, while the percentages with 
respect to the RLM-MAC are 22%, 18%, 28%, 16%, 14% and 2%. Overall, Figures [T] and [3] also suggest 
that in some cases the RLM-MAC will have RMSEs lower than the minimum one of the aLM-EnRML. 

Figure |4] shows the box plots of the data mismatch of the aLM-EnRML (left) and RLM-MAC (right) at 
different iteration steps in one repetition run. The ensemble average data mismatch in both iES decreases 
with the iteration. In this particular case, however, the data mismatch of the RLM-MAC tends to converge 
faster than that of the aLM-EnRML. In addition, its final data mismatch (up to the maximum iteration 
step in the figure) is substantially lower than that of the aLM-EnRML. We also note that, the top outlier 
in each box plot of the aLM-EnRML does not appear to follow the tendency of changes of other data 
mismatch values in the same box plot, which is a phenomenon not observed in the RLM-MAC. A possible 
explanation of this difference might be as follows: In the RLM-MAC, in order for the cost function in 
the MAC problem to decrease with respect to the change of each ensemble member (see the deduction 
in Appendix [Bj, it is natural to use the deviations 

g (m}) - g (rh’j (24) 

to construct the data square root matrix on top of the first order Taylor approximation. In the aLM- 
EnRML, however, AMD is replaced by 

g ( m j) - g (mj) . (25) 

The difference between (l24l) and (l25l) would then depend on the degree of nonlinearity in g. In this 
particular case, because of the strong nonlinearity in the L96 model, the difference between C5l) and (PHI) 
may be relatively large (see Figure [6] later for a numerical investigation), such that the iteration formula 
based on the choice (1251) does not form gradient descent directions for all the ensemble members in the 
aLM-EnRML. 

In accordance to the results in Figure [4l top panels of Figure [5] show the box plots of RMSEs of 
the ensemble members with respect to the initial conditions at different iteration steps in the aLM- 
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EnRML (left) and RLM-MAC (right), and bottom panels indicate the associated ensemble spread^. 
For both the aLM-EnRML and RLM-MAC, their RMSEs drop at the first few iteration steps, bounce 
back subsequently, and then drop again until they converge (although a slight oscillation can also be 
observed around iteration step 85 in the aLM-EnRML). In this particular case, the final RMSEs of the 
RLM-MAC tend to be lower than those of the aLM-EnRML. The ensemble spread of the aLM-EnRML 
monotonically decreases with the iteration, while that of the RLM-MAC reaches its minimum around 
iteration step 10, and then slightly increases afterwards. The final ensemble spreads of both algorithms 
are close to each other, and they are both underestimated in comparison with the corresponding RMSEs. 
One possible reason of this underestimation phenomenon may be because in both smoothers the ensemble 
members tend to converge to regions surrounding certain local minima, therefore the ensemble spreads 
are relatively small, while the RMSEs are determined by the distances between the local minima and 
the truth, and thus are not necessarily as small as the corresponding ensemble spreads. Apart from this, 
we also envision that the effects of finite ensemble size (e.g., sampling errors and correlations among 
ensemble members E3) may also contribute to the underestimation phenomenon. 

Figure H] depicts the normalized differences C^ 1 ^ 2 ^g (m j — g (m®) j at different iteration steps, based 
on the corresponding ensemble members of the RLM-MAC. Please note that in this case, the box plots 
indicate the spatial distribution of the elements of the vector C d 1 ^ 2 ^g (m®) — g (rn® ) j, whose size is 
equal to the number of observations (rather than the ensemble size). In addition, in the box plots, g (in') 


and g (m® ) are computed from the same ensemble of the RLM-MAC at each iteration step, rather than 


one (e.g., g (m 1 )) from the ensemble of the RLM-MAC and another (e.g., g (m® )) from the corresponding 
ensemble of the aLM-EnRML. 

From Figure [6] one can see that the spreads of the normalized differences are relatively large at the 
first few iteration steps. This may be because the initial ensemble spreads are relatively large (see the 
lower right panel of Figure [5]), and the chaos nature of the L96 model tends to further amplify them. 
However, as the number of iterations increases, the ensemble members of the RLM-MAC tend to converge 


7 Suppose that {mjjyJh is an ensemble of models with the ensemble size being N e , and that mj jS is the sth element of 
the vector mj. Let rtf. be the sample variance of the sth elements of all ensemble members, then the ensemble spread is 
defined as cr j /m, where m is the size of the ensemble members mj. 
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toward a certain local minimum, and the ensemble spread decreases rapidly. As a result, g (m 1 ) approach 


g (m*), such that the normalized differences decrease and tend to center around zero. 


A reservoir facies estimation problem 

Here we consider a simple 2D facies model (see Figure [7]) previously studied in [T7] . The reservoir consists 
of 45 x 45 gridblocks, with only oil and water inside. The porosity of the reservoir is 0.1 in each gridblock, 
but the permeability varies spatially: It is 10000 md in the channels, and 500 rnd in the background. In 
the field there are 8 water injectors (II - 18) and 8 producers (PI - P8), and their locations are marked 
as white dots in Figure [3 The injectors are constrained by injection rates (15.9 m 3 /day in each), while 
the producers are constrained by bottom hole pressures (BHP, 399 bar in each). 

The measurements used for history matching are the oil and water production rates in the producers 
and BHP in the injectors. To generate the synthetic measurements, we use the reference facies model in 
Figure [7] as the input to ECLIPSE© [I] and run the simulation for 1900 days. The measurements are 
then taken as the outputs of ECLIPSE© every 190 days, plus certain Gaussian noise (whose variances are 
0.16 m 3 /day for production rates, and 0.07 bar for BHP). To enhance the connectivities of the estimated 
channels, eight statistical measures used in m are also assimilated into the model. These statistics 
measurements include the sample mean values of: (1) area of bodies (a body is defined as neighboring 
connected gridblocks with the same facies type, see [T7]); (2) volume of bodies; (3) area-to-volume ratio 
of bodies; (4) maximal body extension in x-direction; (5) maximal body extension in y-direction; (6) 75th 
percentile for the extensions in the x-direction; (7) number of bodies; and (8) fraction of the total volume of 
the facies with respect to the total volume of the field, where all the mean values are evaluated with respect 
to the initial background ensemble m- Accordingly, in the course of assimilation, the measurement error 
variances associated with these statistics measurements are taken as their sample variances that are also 
evaluated with respect to the initial background ensemble. We have also run both the aLM-EnRML and 
RLM-MAC when the statistics measurements are not included in assimilation. In this case, the final data 
mismatch values are close to what we show in Figure [TU] later, but the corresponding estimated channels 
tend to be less connected (results not shown), consistent with the observation in [17 . 
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We follow [T7] and use the signed distance function in the level set method to characterize the facies. 
When using an iES for history matching, there is no need to augment the model and the corresponding 
state variables as in the EnKF. Therefore for the problem considered here, the initial ensemble only 
contains the signed distance function values at each gridblock, which are converted from the facies models 
generated by SNESIM [35; through the MATLAB level set toolbox (version 1.1, [25]). As in the previous 
example, the ensemble size in the aLM-EnRML is 100, while that in the RLM-MAC is 99. In addition, 
the aforementioned parameter rule 7 * = a 1 trace(S^ (S l d ) T )/N e is also adopted in this example. The 
maximum number of iterations is 20 , and the coefficient /?„ = 2 . 

A few initial models are shown in the first column of Figure [ 8 ] (a, d, g), and the corresponding models 
obtained at the final iteration steps of the aLM-EnRML and RLM-MAC are plotted in the second 
column (b, e, h) and the third one (e, f, i), respectively. As can be seen in these images, despite the clear 
difference among the initial models, all the updated models of the aLM-EnRML and RLM-MAC bear 
certain similarities to each other, and they also, to some extent, resemble the reference model in Figure 
□ 

Figure [9] reports the average scores of the initial ensemble (panel (a)) and the final ensembles of the 
aLM-EnRML (panel (b)) and RLM-MAC (panel (c)) in matching the facies of the reference model. The 
score is calculated as follows: at each gridblock, if a model matches the facies, and it has a score of 1 , 
otherwise it has a score of 0. The final score is then averaged over all the models in the ensemble. Therefore 
a higher score means that more models have the correct estimated facies. From Figure [9] it is seen that, 
compared with the scores in the initial ensemble, it seems that the aLM-EnRML and RLM-MAC tend to 
improve the facies estimation in the channels, but also have deteriorated matching in certain areas (e.g., 
that around the injector 12). In addition, in this particular case, it seems that the RLM-MAC tends to 
perform better than the aLM-EnRML in terms of matching the facies. 

Figure flOl indicates the data mismatch of the aLM-EnRML (left) and RLM-MAC (right) as functions 
of the number of iteration steps. The aLM-EnRML stops at the 10th iteration step due to the stopping 
condition (c) mentioned early (relative change of data mismatch less than 0.01%), while the RLM-MAC 
stops at the maximum iteration step. Despite this difference, in this case the aLM-EnRML and RLM- 
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MAC appear to have comparable performance, in terms of the convergence speed and the final data 
mismatch value. For the aLM-EnRML, its final data mismatch is 2197.71 ± 234.63, where the number 
before ± is the mean, and that after ± is the standard deviations (STD). On the other hand, for the 
RLM-MAC, its final data mismatch is 1864.27 ± 98.9975, slightly lower than that of the aLM-EnRML. 

Following [T7], Figures [TT] and [12] show the performance of history matching oil and water production 
rates in producers P5 and P7, for the initial ensemble and the ensembles of the aLM-EnRML and RLM- 
MAC at their final iteration steps, respectively. In both figures, compared to the results with respect to 
the initial ensemble, those obtained from the final ensembles of the aLM-EnRML and RLM-MAC appear 
to have reduced ensemble spread and improved data mismatch. The aLM-EnRML performs better than 
the RLM-MAC in history matching producer P5, but the situation is reverted in producer P7. In the 
initial ensemble the oil production rates in P5 and the water production rates in P7 tend to be lower than 
the reference values (which may be related to the bias in the facies distribution of the initial ensemble, 
as highlighted in Figure |9 (a) [ ), and this bias trend remains in the final ensembles obtained by the aLM- 
EnRML and RLM-MAC. Similarly, in the initial ensemble the oil production rates in P7 and the water 
production rates in P5 tend to be higher than the reference values. This bias trend is clearly visible in 
the final ensemble of the aLM-EnRML, but seems to be better mitigated in the final ensemble of the 
RLM-MAC. 

The top panels of Figure [T51 show the box plots of RMSEs (in permeability values) of ensemble models 
with respect to the reference model at different iteration steps in the aLM-EnRML (left) and RLM-MAC 
(right), and bottom ones indicate the corresponding spreads in the ensemble. Similar to the situation 
in Figure (5[ oscillations of the RMSEs (e.g., in terms of the medians of the box plots) are observed in 
both iES before the RMSE values enter certain plateaus. In this particular case, the final RMSEs of the 
RLM-MAC tend to be lower than those of the aLM-EnRML. On the other hand, the ensemble spread of 
the aLM-EnRML tends to decrease with the iteration, while that of the RLM-MAC exhibits slight U-turn 
behaviour twice before it enters a plateau. The final ensemble spread of the RLM-MAC is higher than 
that of the aLM-EnRML, and both are underestimated in comparison with the corresponding RMSEs. 
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Figure 1771 shows the normalized differences C^ 1//2 (m'j — g(m*)^ at different iteration steps, 
calculated based on the ensemble models of the RLM-MAC. Compared to Figure [6l it is seen that the 
normalized differences remain relatively small over the iterations, which may be because in this case the 
phase flows are non-turbulent, and they tend not to amplify the differences among the ensemble models 
(in contrast to the turbulence nature in the L96 model of the previous subsection). 


History matching study in the Brugge field case 


Here we compare the history matching performance of the aLM-EnRML and RLM-MAC in the Brugge 
benchmark case study m- The simulation model of the Brugge field case has 60048 (139 x 48 x 9) 
gridblocks, with 44550 of them being active. In the case study there are 20 producers and 10 water 
injectors. The dataset we used is from TNO m and contains 20 years’ production datcjj. The first 
decade’s data are collected at 127 time instants, while the second decade’s at 185 time instants. Following 
we use the production data at 20 out of 127 time instants in the first decade for history matching, 
and the rest of the 20 years’ production data for out-of-sample cross validation, including, for example, 
predicting the oil production rates during the second decade based on the history-matched models. 

In the course of history matching, the producers and injectors are controlled by the liquid rate (LRAT) 
target. The production data consist of the historical oil production rates (WOPRH) and water cuts 
(WWCTH) at 20 producers, and the historical bottom hole pressure (WBHPH) at all 30 wells. There¬ 
fore the total number of production data is 1400. The standard deviations for WOPRH, WWCTH and 
WBHPH are 100 stb/day, 0.05 and 50 psi, respectively. The model variables to be updated include 
porosity (PORO), permeability (PERMX, PERMY, PERMZ) and net-to-gross (NTG) ratio at all ac¬ 
tive gridblocks. Consequently the total number of parameters is 222, 750. In the aLM-EnRML we use 
the 104 model realizations from the initial ensemble provided by TNO E5, while in the RLM-MAC 
we only use 103 of them, after dropping an ensemble member at random (so that the aLM-EnRML 
and RLM-MAC have the same number of forward runs at each iteration). In both iES, the parameter 


8 The second decade’s production data are provided by TNO based on IRIS’ proposed production optimization strategy 


for that period, see, for example, m- 
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rule is 7 * = cP trace(S^(S l d ) T )/N e: with the maximum number of iterations being 20, and the coefficient 

= 0.5. In all experiment runs we use ECLIPSE© [Tj for reservoir simulation. 

Figure [15] indicates the distributions of PERMX on layers 1 and 2. The top images are from a 
realization of the reservoir models in the initial ensemble, the middle ones from the corresponding reservoir 
model updated by the aLM-EnRML, and the bottom ones from the corresponding model obtained by the 
RLM-MAC. Although the final data mismatch values of the aLM-EnRML and RLM-MAC are close to 
each other (see Figure fl6l) . there do not seem to be substantial similarities between the resulting PERMX 
images in Figure 1151 This may be due to the difference of both iES in their ways to approximate the 
gradient of the objective function, hence the difference in the subsequent search directions, as we have 
discussed previously. 

Figure [TC] shows the box plots of data mismatch of the aLM-EnRML (left) and RLM-MAC (right) 
as the function of iteration steps. The aLM-EnRML stops at the 15th iteration step due to the stopping 
condition on the relative change of data mismatch (less than 0.01%), and the RLM-MAC stops at the 
maximum iteration step. The aLM-EnRML and RLM-MAC have close performance again, in terms of 
the convergence speed and the final data mismatch value. In this case, the final data mismatch of the 
aLM-EnRML is 869.90 ± 23.38, while that of the RLM-MAC is 928.96 ± 25.76, about 7% higher in the 
average data mismatch value. 

Figure [IT] depicts the history matching profiles of oil production rates (WOPR) at producers BR- 
P-9 (top row), BR-P-13 (middle row) and BR-P-19 (bottom row) in the first 10 years, with respect to 
the initial ensemble of reservoir models (1st column), the history matched reservoir models of the aLM- 
EnRML using production data at 20 time instants (2nd column), and the counterpart history matched 
reservoir models of the RLM-MAC (3rd column), respectively. In the pictures, the red dots represent the 
historical WOPRH data at the producers, and the blue curves are the forecasts with respect to the initial 
ensemble (1st column), and the history matching profiles of both smoothers (2nd and 3rd columns). 
Consistent with the results in Figure HU the aLM-EnRML and RLM-MAC exhibit close performance in 
matching the WOPRH data. Similar results are also obtained for history matching the WWCTH data at 
BR-P-9, BR-P-13 and BR-P-19, as can be seen in Figure [151 
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Figure Uni also shows the normalized differences C d 1 2 ^g (m 1 ) — g (m* )^ at different iteration steps, 
calculated based on the ensemble models of the RLM-MAC, in the Brugge field case. Similar to Figure 
ns the normalized differences remain relatively small over the iterations, in contrast to the corresponding 
results in the L96 model (see Figure |G|). This may be also because in this case the reservoir fluid dynamics 
is non-turbulent, and tends to be less nonlinear than that in the L96 model. 

Next we present some cross validation results. Figure [20] shows the box plots of the data mismatch 
values in 20 years, with respect to the reservoir models from the initial ensemble (left), the final ensemble 
of the aLM-EnRML using production data at 20 time instants in the first 10 years (middle), and the 
corresponding final ensemble of the RLM-MAC (right). In this case, the data mismatch values with 
respect to the initial ensemble and the final ensembles of the aLM-EnRML and RLM-MAC are 3.48 x 
10 5 ± 2.46 x 10 5 , 4707.93 ± 357.83 and 4540.31 ± 285.34, respectively. Therefore, the reservoir models 
history matched by both iES lead to substantially lower data mismatch than those of the initial ensemble. 
The data mismatch values of the aLM-EnRML and RLM-MAC are close to each other again, with that 
of the RLM-MAC being slightly lower. 

Figure Ell depicts the WOPR profiles at producers BR-P-9 (top row), BR-P-13 (middle row) and 
BR-P-19 (bottom row) in 20 years, with respect to the reservoir models from the initial ensemble (1st 
column), the final ensemble of the aLM-EnRML (2nd column), and the corresponding final ensemble of 
the RLM-MAC (3rd column), respectively. In the pictures, the red dots represent the historical WOPRH 
data at the producers, and the blue curves are the forecasts with respect to the initial ensemble and 
the final ensembles obtained by both smoothers. The pictures in the 2nd and 3rd columns also contain 
some vertical dashed lines, which are used to separate the production periods between the first and 
second decades. In this case, the aLM-EnRML and RLM-MAC again exhibit close performance in cross¬ 
validating the WOPRH data at BR-P-9. However, in the second decade, the aLM-EnRML predicts the 
WOPRH data at BR-P-19 better than the RLM-MAC does, which may be related to the fact that the 
WWCTH data at BR-P-19 is better predicted by the aLM-EnRML (see Figure l22l) . The situation is 
reverted for the prediction of the WOPRH data at BR-P-13 in the second decade, and the RLM-MAC 
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outperforms the aLM-EnRML instead. This may be because the WWCTH data at BR-P-13 is better 
predicted by the RLM-MAC instead, as can be seen in Figure 1221 


Conclusion 

In this work we showed that an iteration formula similar to those used in the ensemble smoother with 
multiple data assimilation (ES-MDA) and the approximate Levenberg-Marquardt ensemble randomized 
maximum likelihood (aLM-EnRML) method can be derived from the regularized Levenberg-Marquardt 
(RLM) algorithm in inverse problems theory, when history matching is recast as a minimum-average-cost 
(MAC) problem. 

Specifically, to reduce the computational cost in solving the MAC problem, a simple strategy is to 
approximate the nonlinear simulator, for all ensemble models, through a first order Taylor expansion 
around a single common point, e.g., the ensemble mean at the previous iteration step as considered in 
the current study. In addition, to avoid directly evaluating the resulting (common) Jacobian matrix, one 
may follow the convention in the ensemble Kalman filter (EnKF) to approximate the product between 
the Jacobian matrix and the model square root matrix by the corresponding data square root matrix. In 
this way, the resulting iteration formula is similar to those in the ES-MDA and aLM-EnRML, but differs 
in the way of constructing the data square root matrix: In the ES-MDA and aLM-EnRML, the simulated 
data center around their sample mean, while in the RLM-MAC, they center around the simulated data 
of the mean of the ensemble models. In general, it is expected that one may solve the MAC problem in 
a more accurate way - but possibly with a higher computational cost - by evaluating the gradients at 
multiple common points. 

To examine the effect of different ways in constructing the data square root matrix, we compare the 
performance of the aLM-EnRML and RLM-MAC in three numerical examples. In the first example we 
consider a strongly nonlinear system (the Lorenz 96 model). In this particular case the aLM-EnRML and 
RLM-MAC exhibit more diverged behaviour and the RLM-MAC tends to perform better than the aLM- 
EnRML. In the second and third examples we apply the aLM-EnRML and RLM-MAC to a reservoir 
facies estimation problem and the history matching problem in the Brugge field case, respectively. In 
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both cases it appears that the aLM-EnRML and RLM-MAC exhibit close performance due to weaker 
nonlinearity in the reservoir fluid dynamics. 
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Nomenclature 

d = data 

d° = observation data 

f = function 

g = reservoir simulator 

m = dimension of reservoir model 

m = reservoir model 

m c = common reservoir model in the minimum-average-cost problem 

m = ensemble-average reservoir model 

N e = ensemble size 

p = number of observation data 

Pk = probability density function (pdf) at the fcth iteration step 
r = scalar coefficient larger than 1 

t = unit-less sudo-time variable in the Lorenz 96 model 
v = generic random vector 

v = estimation of v 

v tr = true realization value of the random vector v 

x = state variable in the Lorenz 96 model 

y = random vector that represents observations 

A = generic matrix 

C d = observation error covariance matrix 

C m = model error covariance matrix as in the EnKF 
C m = modified model error covariance matrix in the RLM-MAC 
F = unit-less driving force in the Lorenz 96 model 
1 ~ = functional 

Q = inverse functional of F 

G = Jacobian of reservoir simulator 

Q c = composite functional 

/ = number of iterations in ES-MDA 
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I p = identity matrix of size p 

M = ensemble of reservoir models 

O = objective function 

6 = expectation of O 

S d = data square root matrix as in the EnKF 
Srf = modified data square root matrix in the RLM-MAC 
S m = model square root matrix 


Greek symbols 

a = adaptive coefficient in the RLM-MAC that determines 7 

6 = Dirac delta function 

7 = adaptive coefficient in iterative ensemble smoothers 

er 2 = sample variances of the ensemble of reservoir models 
p = tapering factor of 7 

Am = model change 


Subscript 


a = analysis 

b = background 

d = data 

e = ensemble 

j = index of ensemble member 

k = iteration index 

to = model 

s = index of sample variances of the ensemble of reservoir models 
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Superscript 

b = background 

i = iteration index 

o = observation 

tr = truth 



28 


Xiaodong Luo et al. 


Acknowledgements We would like to thank Drs. Yan Chen (Total), Randi Valestrand (IRIS) and three anonymous 
reviewers for their valuable helps, comments and suggestions, which have significantly improved the work. We also acknowl¬ 
edge the IRIS/CIPR cooperative research project “Integrated Workflow and Realistic Geology” which is funded by industry 
partners ConocoPhillips, Eni, Petrobras, Statoil, and Total, as well as the Research Council of Norway (PETROMAKS) for 
financial support, and thank Schlumberger for ECLIPSE academic licenses and the developers of SGeSM and MATLAB 
level set toolbox (Ian M. Mitchell) for providing their softwares. 


A Characterizing the RLM-MAC from the point of view of an expectation-maximization 
algorithm 


Following a recent working paper m, the RLM-MAC algorithms can also be described in analogy to the expectation- 
maximization (EM) algorithm [5], which is an iterative algorithm that consists of an expectation step (E-step) and an 
maximization step (M-step) at each iteration. At the E-step, one constructs an expectation-type cost function that depends 
on certain unknown quantities-of-interest (QOI) and their prior knowledge, while at the subsequent M-step, one solves a 
maximization problem to obtain an updated estimation of the QOI. The updated QOI are then passed to the next iteration 
as the new prior knowledge of the corresponding E-step, and so on. 

The cost function in the RLM-MAC can be considered as a Monte-Carlo approximation to the following expectation-type 
cost function at the E-step of the kth iteration: 

6 k = f O k (mjj, mj,y) p k (m^, , y|y) dmj dy , (A.l) 

where is the prior knowledge of the QOI, mjj the QOI to be updated, y the real observation, and y a perturbation of 
y; Ok{ mjj, m^,y) is a chosen function (e.g., as the quadratic function in the RLM-MAC) and may depend on m^, mjj and 
y, and Ok is the expectation of Ofc(m^,m^,y) with respect to the (conditional) joint pdf m^, y|y). In the context 

of data assimilation, the updated QOI mjj often depend on both and y, so that one has mjj = J~k (m^,y), where J~k 
is an unknown deterministic map (or functional), and is determined at the subsequent “M-step” under a certain criterion. 
In the RLM-MAC, Tk are determined by solving the following functional optimization problem 


argmaxOfc . 
Fk 


(A.2) 


Note that J~k may also be chosen from other points of views. For instance, they may be constructed in a way such that 
the conditional pdf pfcfmjjly) is consistent with Bayes’ rule, in the sense of satisfying 


p.( m ;iy)- p(y ™ (mi) . 

p(y) 


(A.3) 


On the other hand, note that = mj Through recursion, one has =T k (^F k —1 ( m j 1 • y) . y j = 

■. Let mj = JFC ( m®, y), where is the composition of the functionals Ti , then by the definition of pdf of transformed 
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random variables we have the relation between po( m 5 ?y|y) and Pfc(m^|y) given by 

d 


Pfc( m Jy) = 




/ 


Po(m£,y|y) drn^dy 


^(mg,y)<m«S 


8 


dm% 


J Po(m£)p(y|y) drn^dy , 


(A.4) 


where p(yjy) is the pdf of the perturbed observation y, given the real observation y. Note that in the second line of Eq. 
JA.4L we have assumed that the initial prior m® is independent of observation. Combining Eqs. dA.41) and (lA.3l) , we require 
the following identity 

d 




( , On x j 0j~ P(y| m S)Po(mJ) 

J Po(m b )p(y|y) dm b dy = -- 


(A.5) 


(mg,y)<mfc 


to hold in order for p^{ m*|y) to satisfy Bayes’ rule. 

Eq. (IA.5I) suggests that the choice of Ti depends on both the prior pdf po(•) and the likelihood p(y|»). As a side remark, 
we also note that Eq. itAAl only imposes a constraint on the compositional functional , but does not specify how J~i 
should be chosen at individual iteration steps. This freedom allows one to incorporate additional criteria (e.g., in terms of 
solving optimization problems in iron at some of the iteration steps, subject to the condition in Eq. JAA]). 

In linear systems with Gaussian model and observation errors, [6] provides proof and examples on how to construct the 
functionals to make the resulting p^( mjj|y) consistent with Bayes’ rule. In general situations with nonlinearity and/or 
non-Gaussianity, the constructions of become more complicated and are beyond the scope of the current study. 


B Derivation of the iteration formula in RLM-MAC 


By inserting Eq. cei into G2), one has the following approximate cost function 


Ji+1 =■$- E [ ( d ° - s K) - g = B +1 - m 0) T c 4 1 ( d ? - g K) - g = B +1 - m 0) (B - 1} 

e 3 =1 


m’ + 1 -m'| ( Cl, 


B +1 -”5 


A necessary condition for mV* -1 (j = 1,2, , N e ) being the solution of the MAC problem is thus dJ t+1 /d m“ +1 = 0, 

which leads to the following equation: 


-2 ( G T G d 1 ( d r - 8 K) - G* (m; +1 - m*)) + 2 7 ‘ (cjj (mf 1 - m} j = 0. 


(B.2) 


With some algebraic operations, the above equation becomes 


(G*) T C d 1 G ! c + Y (cQ ’l (m} +1 - mj) = (G*) T C d 1 [dj - g (m l c ) - G* (mj - m*)] . (B.3) 


By the Sherman-Morrison-Woodbury formula m, one has 


(G”) T C d - 1 G‘+ 7 i (C* 


(G*) C d 1 = (Gc) GjCj„(G') +Yc d 


Inserting the above identity into Eq. IIB.3I ) and letting i one obtains Eq. 11191 . 
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List of Figures 


1 Histograms of (a) data mismatch and (b) RMSEs of the aLM-EnRML and RLM-MAC obtained from 100 
repetitions of the experiment in the L96 model. For visualization, the horizontal axes of Figure [i(a)| are in the 


logarithmic scale. In Figure [1(b) | the RMSEs are calculated based on the differences between the trajectory 
generated by the true initial state (in 10 days) and those generated by the estimated ones. 

2 Pie charts of the data mismatch of the aLM-EnRML (left) and RLM-MAC (right) with respect to the results 

in Fig. [I] . |3. 

3 Pie charts of the RMSEs of the aLM-EnRML (left) and RLM-MAC (right) with respect to the results in Fig. ^ [3; 

4 Box plots of the data mismatch of the aLM-EnRML (left) and the RLM-MAC (right) at different iteration 

steps in one run in the L96 model. The maximum number of iterations for the RLM-MAC is 100, while 
that for the aLM-EnRML is increased to 200 and the aLM-EnRML converges after 120 iteration steps. For 
visualization, the vertical axes are in the logarithmic scale, and the results are plotted at iteration steps 0, 5, 
and every 5 iteration steps thereafter, along the horizontal axes. In each box plot, the horizontal line (in red) 
inside the box denotes the median; the top and bottom of the box represent the 75th and 25 th percentiles, 
respectively; the whiskers indicate the ranges beyond which the data are considered outliers, and the whiskers’ 
positions are determined using the default MATLAB© setting (please see the corresponding documentation 
in MATLAB© R2012a), while the outliers themselves are plotted individually as plus signs (in red). 


Top panels: Box plots of RMSEs of the ensemble members with respect to the initial conditions at different 
iteration steps in the aLM-EnRML (left) and RLM-MAC (right). Bottom panels: Ensemble spreads of the 
aLM-EnRML (left) and RLM-MAC (right) at different iteration steps of the L96 model. 


6 Box plots of the normalized differences C d 
model. 


1/2 


(m 1 ) — g J at different iteration steps in the L96 
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7 Reference facies model. The white dots indicate the locations of injectors (II - 18) and producers (PI - P8). |42 

8 A few models (members 1 - 3) in the initial ensemble (1st column) and the ensembles of the aLM-EnRML 
(2nd column) and RLM-MAC (3rd column) at the final iteration steps of the facies estimation problem. . . 

9 Average facies matching scores of (a) the initial ensemble, (b) the final ensemble of the aLM-EnRML and (c) 

the final ensemble of the RLM-MAC. A higher score means that more models have the correct facies. |44 


10 Box plots of data mismatch at different iteration steps of the facies estimation problem. Left: aLM-EnRML; 

Right: RLM-MAC. |4 
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11 History matching of oil production rates at P5 (top row) and P7 (bottom row) of the facies estimation 
problem, using the initial ensemble (1st column) and the ensembles of the aLM-EnRML (2nd column) and 
RLM-MAC (3rd column) at the final iteration steps. In each image, the red dots represent the synthetic 
measurements and the blue curves are the forecasts with respect to the initial ensemble (1st column) and 
history matching profiles (2nd and 3rd columns). For comparison, the production data of the reference model 
are also shown in orange. 


12 As in Figure ITTl but for history matching of water production rates in the facies estimation problem. 

13 Top panels: Box plots of RMSEs of the ensemble models—with respect to the reference model of the facies 
estimation problem-at different iteration steps in the aLM-EnRML (left) and RLM-MAC (right). Bottom 
panels: Ensemble spreads of the aLM-EnRML (left) and RLM-MAC (right) at different iteration steps. . . . 
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Box plots of the normalized differences 
estimation problem. 



at different iteration steps of the facies 


15 Distributions of PERMX on layers 1 and 2 from a reservoir model of the initial ensemble (a,b), the corre¬ 
sponding reservoir model updated by the aLM-EnRML (c,d), and the corresponding reservoir model updated 
by the RLM-MAC (e,f). The dots in the figures indicate the locations of the producers and injectors in the 
Brugge field. 


y 


16 Box plots of data mismatch of the aLM-EnRML (left) and the RLM-MAC (right) at different iteration steps 
in the Brugge field case, using production data at 20 out of 127 time instants in the first 10 years. 


51 


17 History matching of oil production rates (WOPR) at BR-P-9 (top), BR-P-13 (middle) and BR-P-19 (bottom) 

in the Brugge field case, using the initial ensemble (1st column) and the ensembles of the aLM-EnRML (2nd 
column) and RLM-MAC (3rd column) at the final iteration steps. In each image, the red dots represent 
the historical WOPRH data and the blue curves are the forecasts with respect to the initial ensemble (1st 
column), and history matching profiles (2nd and 3rd columns). 

18 As in Figure [TT] but for water cut (WWCT) at BR-P-9, BR-P-13 and BR-P-19 in the Brugge field case. . . 

19 Box plots of the normalized differences C d 1//2 (rind) — g J at different iteration steps in the Brugge 

field case. 


y 

3 


20 Box plots of data mismatch with respect to the initial ensemble (left), the aLM-EnRML (middle) and the 
RLM-MAC (right) in the Brugge field case. Here, all the production data in the first 20 years are used to cross- 
validate the history-matched reservoir models obtained by aLM-EnRML and RLM-MAC using production 
data at 20 out of 127 time instants in the first 10 years. 
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21 Cross validation of oil production rates (WOPR) at BR-P-9 (top), BR-P-13 (middle) and BR-P-19 (bottom) 
in the Brugge field case, using the initial ensemble (1st column) and the ensembles of the aLM-EnRML (2nd 
column) and RLM-MAC (3rd column) at the final iteration steps. In each image, the red dots represent 
the historical WOPRH data and the blue curves are the forecasts with respect to the initial ensemble (1st 
column) and the final ensembles obtained by both iES (2nd and 3rd columns). The vertical dashed lines in 
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the figures of the 2nd and 3rd columns separate the first and second decades. 

22 As in Figure fT71 but for water cut (WWCT) at BR-P-9, BR-P-13 and BR-P-19 in the Brugge field case. . 
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Data mismatch (aLM-EnRML) Data mismatch (RLM-MAC) RMSE (aLM-EnRML) RMSE (RLM-MAC) 


(a) Histogram of average data mismatch (b) Histogram of average RMSEs 

Fig. 1 Histograms of (a) data mismatch and (b) RMSEs of the aLM-EnRML and RLM-MAC obtained from 100 repetitions 
of the experiment in the L96 model. For visualization, the horizontal axes of Figure [l (a) | are in the logarithmic scale. In 
Figure [l (b) | the RMSEs are calculated based on the differences between the trajectory generated by the true initial state 
(in 10 days) and those generated by the estimated ones. 
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(a) aLM-EnRML (b) RLM-MAC 

Fig. 2 Pie charts of the data mismatch of the aLM-EnRML (left) and RLM-MAC (right) with respect to the results in 

Fig.m 
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(a) aLM-EnRML 


(b) RLM-MAC 


Fig. 3 Pie charts of the RMSEs of the aLM-EnRML (left) and RLM-MAC (right) with respect to the results in Fig. \T\ 
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Fig. 4 Box plots of the data mismatch of the aLM-EnRML (left) and the RLM-MAC (right) at different iteration steps in 
one run in the L96 model. The maximum number of iterations for the RLM-MAC is 100, while that for the aLM-EnRML 
is increased to 200 and the aLM-EnRML converges after 120 iteration steps. For visualization, the vertical axes are in 
the logarithmic scale, and the results are plotted at iteration steps 0, 5, and every 5 iteration steps thereafter, along the 
horizontal axes. In each box plot, the horizontal line (in red) inside the box denotes the median; the top and bottom 
of the box represent the 75th and 25 th percentiles, respectively; the whiskers indicate the ranges beyond which the data 
are considered outliers, and the whiskers’ positions are determined using the default MATLAB© setting (please see the 
corresponding documentation in MATLAB© R2012a), while the outliers themselves are plotted individually as plus signs 
(in red). 
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Fig. 5 Top panels: Box plots of RMSEs of the ensemble members with respect to the initial conditions at different iteration 
steps in the aLM-EnRML (left) and RLM-MAC (right). Bottom panels: Ensemble spreads of the aLM-EnRML (left) and 
RLM-MAC (right) at different iteration steps of the L96 model. 

















Title Suppressed Due to Excessive Length 


41 


Fig. 6 


-100M 


-200 t 


-300 


+■ 4- 4- 


10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 

Iteration number (RLM-MAC) 


Box plots of the normalized differences 


r ~ 1/2 



at different iteration steps in the L96 model. 






42 


Xiaodong Luo et al. 


Fig. 7 


Reference facies model 



Reference facies model. The white dots indicate the locations of injectors (II - 18) and producers (PI - P8). 
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(g) Initial ensemble 


(h) aLM-EnRML 


(i) RLM-MAC 


Fig. 8 A few models (members 1 - 3) in the initial ensemble (1st column) and the ensembles of the aLM-EnRML (2nd 


column) and RLM-MAC (3rd column) at the final iteration steps of the facies estimation problem. 
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(a) Initial ensemble (b) aLM-EnRML (c) RLM-MAC 

Fig. 9 Average facies matching scores of (a) the initial ensemble, (b) the final ensemble of the aLM-EnRML and (c) the 
final ensemble of the RLM-MAC. A higher score means that more models have the correct facies. 
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Fig. 10 Box plots of data mismatch at different iteration steps of the facies estimation problem. Left: aLM-EnRML; Right: 
RLM-MAC. 
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(d) Initial ensemble (e) aLM-EnRML (f) RLM-MAC 

Fig. 11 History matching of oil production rates at P5 (top row) and P7 (bottom row) of the facies estimation problem, 
using the initial ensemble (1st column) and the ensembles of the aLM-EnRML (2nd column) and RLM-MAC (3rd column) 
at the final iteration steps. In each image, the red dots represent the synthetic measurements and the blue curves are 
the forecasts with respect to the initial ensemble (1st column) and history matching profiles (2nd and 3rd columns). For 


comparison, the production data of the reference model are also shown in orange. 
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Fig. 12 As in Figure [TT1 but for history matching of water production rates in the facies estimation problem. 
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Fig. 13 Top panels: Box plots of RMSEs of the ensemble models—with respect to the reference model of the facies 
estimation problem-at different iteration steps in the aLM-EnRML (left) and RLM-MAC (right). Bottom panels: Ensemble 
spreads of the aLM-EnRML (left) and RLM-MAC (right) at different iteration steps. 
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Fig. 14 


Box plots of the normalized differences 



at different iteration steps of the facies 


estimation problem. 
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Fig. 15 Distributions of PERMX on layers 1 and 2 from a reservoir model of the initial ensemble (a,b), the corresponding 
reservoir model updated by the aLM-EnRML (c,d), and the corresponding reservoir model updated by the RLM-MAC (e,f). 
The dots in the figures indicate the locations of the producers and injectors in the Brugge field. 






































Title Suppressed Due to Excessive Length 


51 



0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 

Iteration number (aLM-EnRML) 


2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 


(a) aLM-EnRML 


Iteration number (RLM-MAC) 

(b) RLM-MAC 


Fig. 16 Box plots of data mismatch of the aLM-EnRML (left) and the RLM-MAC (right) at different iteration steps in 
the Brugge field case, using production data at 20 out of 127 time instants in the first 10 years. 
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Fig. 17 History matching of oil production rates (WOPR) at BR-P-9 (top), BR-P-13 (middle) and BR-P-19 (bottom) 
in the Brugge field case, using the initial ensemble (1st column) and the ensembles of the aLM-EnRML (2nd column) and 
RLM-MAC (3rd column) at the final iteration steps. In each image, the red dots represent the historical WOPRH data and 
the blue curves are the forecasts with respect to the initial ensemble (1st column), and history matching profiles (2nd and 
3rd columns). 
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Fig. 18 As in Figure [TT] but for water cut (WWCT) at BR-P-9, BR-P-13 and BR-P-19 in the Brugge field case. 
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Fig. 19 Box plots of the normalized differences 


at different iteration steps in the Brugge field 
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Fig. 20 Box plots of data mismatch with respect to the initial ensemble (left), the aLM-EnRML (middle) and the RLM- 
MAC (right) in the Brugge field case. Here, all the production data in the first 20 years are used to cross-validate the 
history-matched reservoir models obtained by aLM-EnRML and RLM-MAC using production data at 20 out of 127 time 
instants in the first 10 years. 
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Fig. 21 Cross validation of oil production rates (WOPR) at BR-P-9 (top), BR-P-13 (middle) and BR-P-19 (bottom) in 
the Brugge field case, using the initial ensemble (1st column) and the ensembles of the aLM-EnRML (2nd column) and 
RLM-MAC (3rd column) at the final iteration steps. In each image, the red dots represent the historical WOPRH data 
and the blue curves are the forecasts with respect to the initial ensemble (1st column) and the final ensembles obtained by 
both iES (2nd and 3rd columns). The vertical dashed lines in the figures of the 2nd and 3rd columns separate the first and 


second decades. 
































































Water cut Water cut Water cut 


Title Suppressed Due to Excessive Length 


57 





(g) Initial ensemble 


(h) aLm-EnRML 


(i) RLM-MAC 


Fig. 22 As in Figure [TT] but for water cut (WWCT) at BR-P-9, BR-P-13 and BR-P-19 in the Brugge field case. 
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